home *** CD-ROM | disk | FTP | other *** search
/ Fritz: All Fritz / All Fritz.zip / All Fritz / FILES / PROGMISC / PCSSP.LZH / PC-SSP.ZIP / MATINSL.ZIP / GELG.FOR < prev    next >
Text File  |  1985-11-29  |  5KB  |  166 lines
  1. C
  2. C     ..................................................................
  3. C
  4. C        SUBROUTINE GELG
  5. C
  6. C        PURPOSE
  7. C           TO SOLVE A GENERAL SYSTEM OF SIMULTANEOUS LINEAR EQUATIONS.
  8. C
  9. C        USAGE
  10. C           CALL GELG(R,A,M,N,EPS,IER)
  11. C
  12. C        DESCRIPTION OF PARAMETERS
  13. C           R      - THE M BY N MATRIX OF RIGHT HAND SIDES.  (DESTROYED)
  14. C                    ON RETURN R CONTAINS THE SOLUTION OF THE EQUATIONS.
  15. C           A      - THE M BY M COEFFICIENT MATRIX.  (DESTROYED)
  16. C           M      - THE NUMBER OF EQUATIONS IN THE SYSTEM.
  17. C           N      - THE NUMBER OF RIGHT HAND SIDE VECTORS.
  18. C           EPS    - AN INPUT CONSTANT WHICH IS USED AS RELATIVE
  19. C                    TOLERANCE FOR TEST ON LOSS OF SIGNIFICANCE.
  20. C           IER    - RESULTING ERROR PARAMETER CODED AS FOLLOWS
  21. C                    IER=0  - NO ERROR,
  22. C                    IER=-1 - NO RESULT BECAUSE OF M LESS THAN 1 OR
  23. C                             PIVOT ELEMENT AT ANY ELIMINATION STEP
  24. C                             EQUAL TO 0,
  25. C                    IER=K  - WARNING DUE TO POSSIBLE LOSS OF SIGNIFI-
  26. C                             CANCE INDICATED AT ELIMINATION STEP K+1,
  27. C                             WHERE PIVOT ELEMENT WAS LESS THAN OR
  28. C                             EQUAL TO THE INTERNAL TOLERANCE EPS TIMES
  29. C                             ABSOLUTELY GREATEST ELEMENT OF MATRIX A.
  30. C
  31. C        REMARKS
  32. C           INPUT MATRICES R AND A ARE ASSUMED TO BE STORED COLUMNWISE
  33. C           IN M*N RESP. M*M SUCCESSIVE STORAGE LOCATIONS. ON RETURN
  34. C           SOLUTION MATRIX R IS STORED COLUMNWISE TOO.
  35. C           THE PROCEDURE GIVES RESULTS IF THE NUMBER OF EQUATIONS M IS
  36. C           GREATER THAN 0 AND PIVOT ELEMENTS AT ALL ELIMINATION STEPS
  37. C           ARE DIFFERENT FROM 0. HOWEVER WARNING IER=K - IF GIVEN -
  38. C           INDICATES POSSIBLE LOSS OF SIGNIFICANCE. IN CASE OF A WELL
  39. C           SCALED MATRIX A AND APPROPRIATE TOLERANCE EPS, IER=K MAY BE
  40. C           INTERPRETED THAT MATRIX A HAS THE RANK K. NO WARNING IS
  41. C           GIVEN IN CASE M=1.
  42. C
  43. C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
  44. C           NONE
  45. C
  46. C        METHOD
  47. C           SOLUTION IS DONE BY MEANS OF GAUSS-ELIMINATION WITH
  48. C           COMPLETE PIVOTING.
  49. C
  50. C     ..................................................................
  51. C
  52.       SUBROUTINE GELG(R,A,M,N,EPS,IER)
  53. C
  54. C
  55.       DIMENSION A(1),R(1)
  56.       IF(M)23,23,1
  57. C
  58. C     SEARCH FOR GREATEST ELEMENT IN MATRIX A
  59.     1 IER=0
  60.       PIV=0.
  61.       MM=M*M
  62.       NM=N*M
  63.       DO 3 L=1,MM
  64.       TB=ABS(A(L))
  65.       IF(TB-PIV)3,3,2
  66.     2 PIV=TB
  67.       I=L
  68.     3 CONTINUE
  69.       TOL=EPS*PIV
  70. C     A(I) IS PIVOT ELEMENT. PIV CONTAINS THE ABSOLUTE VALUE OF A(I).
  71. C
  72. C
  73. C     START ELIMINATION LOOP
  74.       LST=1
  75.       DO 17 K=1,M
  76. C
  77. C     TEST ON SINGULARITY
  78.       IF(PIV)23,23,4
  79.     4 IF(IER)7,5,7
  80.     5 IF(PIV-TOL)6,6,7
  81.     6 IER=K-1
  82.     7 PIVI=1./A(I)
  83.       J=(I-1)/M
  84.       I=I-J*M-K
  85.       J=J+1-K
  86. C     I+K IS ROW-INDEX, J+K COLUMN-INDEX OF PIVOT ELEMENT
  87. C
  88. C     PIVOT ROW REDUCTION AND ROW INTERCHANGE IN RIGHT HAND SIDE R
  89.       DO 8 L=K,NM,M
  90.       LL=L+I
  91.       TB=PIVI*R(LL)
  92.       R(LL)=R(L)
  93.     8 R(L)=TB
  94. C
  95. C     IS ELIMINATION TERMINATED
  96.       IF(K-M)9,18,18
  97. C
  98. C     COLUMN INTERCHANGE IN MATRIX A
  99.     9 LEND=LST+M-K
  100.       IF(J)12,12,10
  101.    10 II=J*M
  102.       DO 11 L=LST,LEND
  103.       TB=A(L)
  104.       LL=L+II
  105.       A(L)=A(LL)
  106.    11 A(LL)=TB
  107. C
  108. C     ROW INTERCHANGE AND PIVOT ROW REDUCTION IN MATRIX A
  109.    12 DO 13 L=LST,MM,M
  110.       LL=L+I
  111.       TB=PIVI*A(LL)
  112.       A(LL)=A(L)
  113.    13 A(L)=TB
  114. C
  115. C     SAVE COLUMN INTERCHANGE INFORMATION
  116.       A(LST)=J
  117. C
  118. C     ELEMENT REDUCTION AND NEXT PIVOT SEARCH
  119.       PIV=0.
  120.       LST=LST+1
  121.       J=0
  122.       DO 16 II=LST,LEND
  123.       PIVI=-A(II)
  124.       IST=II+M
  125.       J=J+1
  126.       DO 15 L=IST,MM,M
  127.       LL=L-J
  128.       A(L)=A(L)+PIVI*A(LL)
  129.       TB=ABS(A(L))
  130.       IF(TB-PIV)15,15,14
  131.    14 PIV=TB
  132.       I=L
  133.    15 CONTINUE
  134.       DO 16 L=K,NM,M
  135.       LL=L+J
  136.    16 R(LL)=R(LL)+PIVI*R(L)
  137.    17 LST=LST+M
  138. C     END OF ELIMINATION LOOP
  139. C
  140. C
  141. C     BACK SUBSTITUTION AND BACK INTERCHANGE
  142.    18 IF(M-1)23,22,19
  143.    19 IST=MM+M
  144.       LST=M+1
  145.       DO 21 I=2,M
  146.       II=LST-I
  147.       IST=IST-LST
  148.       L=IST-M
  149.       L=A(L)+.5
  150.       DO 21 J=II,NM,M
  151.       TB=R(J)
  152.       LL=J
  153.       DO 20 K=IST,MM,M
  154.       LL=LL+1
  155.    20 TB=TB-A(K)*R(LL)
  156.       K=J+L
  157.       R(J)=R(K)
  158.    21 R(K)=TB
  159.    22 RETURN
  160. C
  161. C
  162. C     ERROR RETURN
  163.    23 IER=-1
  164.       RETURN
  165.       END
  166.